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ABSTRACT 

We numerically calculate the growth and saturation of the Rayleigh- Taylor instability caused by 
the deceleration of relativistic outflows with Lorentz factor V = 10, 30, and 100. The instability 
generates turbulence whose scale exhibits strong dependence on Lorentz factor, as only modes within 
the causality scale A9 ~ 1/r can grow. We develop a simple diagnostic to measure the fraction 
of energy in turbulent eddies and use it to estimate magnetic field amplification by the instability. 
We estimate a magnetic energy fraction eb ~ 10 -2 due to Rayleigh- Taylor turbulence in a shock- 
heated region behind the forward shock. The instability completely disrupts the contact discontinuity 
between the ejecta and the swept up circumburst medium. The reverse shock is stable, but is impacted 
by the Rayleigh- Taylor instability, which strengthens the reverse shock and pushes it away from the 
forward shock. The forward shock front is unaffected by the instability, but Rayleigh- Taylor fingers 
can penetrate about 10% of the way into the energetic region behind the shock during the two-shock 
phase of the explosion. We calculate afterglow emission from the explosion and find the reverse shock 
emission to be significantly altered by the instability. The reverse shock emission peaks at a later time 
but is still distinguishable from the forward shock. These calculations are performed using a novel 
numerical technique that includes a moving computational grid. The moving grid is essential as it 
maintains contact discontinuities to high precision and can easily evolve flows with Lorentz factors 
upwards of 300. 

Subject headings: hydrodynamics - methods: numerical - relativity 



1. INTRODUCTION 

In the fireball model (jPaczvnski I 119861 : iGoodmanl 
[Mi), gamma ray bursts are thought to be generated 
by a hot explosion which expands and compresses itself 
into a thin ultrarelativistic shell. Internal collisions of 
such shells are thought to generate the prompt burst of 
gamma rays, after which afterglow emission is produced 
during the further expansion of the shell. Eventually 
this shell transfers its energy into a relativistic blastwave 
propagating into the circumburst medium. During the 
transfer process there is an interesting phase of the evo- 
lution, during which the shell is unstable to the Rayleigh- 
Taylor instability. 

The instability requires a density gradient, for example 
at the interface between two fluids of different density. In 
this case, the two fluids are the material contained in the 
original explosion (the ejecta) and the ambient gas swept 
up by the explosion (the circumburst medium). Once the 
explosion has begun to sweep up a non-negligible amount 
of mass, two shocks form. One of the shocks drives its 
way forward into the ambient gas, and the other back into 
the ejecta. Between the two shocks resides the contact 
discontinuity, marking the separation between ambient 
gas and ejecta. It is this contact discontinuity which is 
unstable. 

The instability should have some observable signa- 
tures. This is primarily due to its impact on the reverse 
shock, which might be observable in the afterglow ra- 
diation from the burst, or perhaps even in the prompt 
emission. The dynamics of the Rayleigh- Taylor insta- 
bility should be taken into account in order to accu- 
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rately predict the magnitude and frequency of the so- 
called "optical flash" associated with the reverse shock. 
Additionally, if the instability generates enough shear to 
cascade into Kelvin-Hclmholtz turbulence, this could am- 
plify magnetic fields which facilitate synchrotron emis- 
sion. Finally, it is unknown whether Rayleigh- Taylor 
fingers mig ht propagate all the way up to the for- 
ward shock (jLevinson Il2009al ). possibly altering the fun- 
damental dynamics of the shock, which arc generally 
thou ght to be governed by the Blandford-McKce solu- 
tion ( Blandford fe McKee I Il976f ) in the ultrarelativistic 
limit. 

These considerations motivate the calculation of the 
growth of the Rayleigh- Taylor instability in the physical 
context of a relativi stic fireball. The nonrelativistic case 
was first treated bv I Chevalier et aTl ()1992j ). who analyt- 
ically calculated the linear growth rate of the instabil- 
ity, and also explored numerically its l i near g rowth and 
nonlinear saturation. Uun fe Norman I (|1996h later per- 
formed a two-dimensional magnetohydrodynamics calcu- 
lation which demonstrated how magnetic fields tend to 
align themselves along Rayleigh- Taylor fingers. 

The relativistic case is relevant for gamma ray bursts, 
which have Lorentz factors T > 100, and relativistic ef- 
fects should qualitatively change the physics of the in- 
stability. First, the growt h rate is modified by time 
dilation (jWang et al 1 120021 ). Secondly, causality argu- 
ments dictate which modes grow and saturate during 
different phases of evolution, analogous to the growth 
of structure in cosmology. Additionally, any turbu- 
lence being generated is relativistic turbulence, which 
may hav e distinct properties from nonrelativistic tur- 
bulence (jZrake fc MacFadyen lf2012f ). Finally, while the 
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Fig. 1. — The structure of the turbulence strongly depends on relativistic effects. When the energy per mass T is varied (the characteristic 
Lorentz factor of the initial shell), this affects the scale and prominence of Rayleigh- Taylor turbulence. This is a snapshot at t = for 
r = 10, 30, and 100. On the lower part of the figure, we plot logarithm of density, while on the upper part we plot a passive scalar which 
tracks the growth of the instability. The dominant mode is of order the relativistic causality scale A8 ~ 1/r. 



Sedov- Taylor solution is a fast attractor, the relativistic 
Blan dford-McKcc solution is a slow attractor ()Gruzinov I 
120001 ). Any deviation from Blandford-McKee due to 
Rayleigh- Taylor fingers might persist until the shock be- 
comes nonrelativistic. 

To exten d the nonrelativis tic results into the relativis- 
tic regime, iLevinsonl (|2009bD performed a stability anal- 
ysis o n the two-shock solution ([Nakamura fc Shigevamal 
I2006D and found linear growth rates which could poten- 
tially be large enough to impact the forward shock. 

In this work, we evolve a relativistic fireball numeri- 
cally in 2D, to calculate its dynamics, and to explore the 
consequences of this instability to observations of gamma 
ray bursts and their afterglows. We study in detail how 
the qualitative features of turbulence are changed by rel- 
ativistic effects (FigureQ}, and we describe generally how 
the dynamics of the two-shock system are affected. To 
improve code performance, we use a moving numerical 
mesh, which follows the radial motion of the bulk flow 
and reduces the advection errors associated with this mo- 
tion. Doing this results in a vastly improved calculation 
of the shock dynamics, and an increased sensitivity to 
the subtle effects that must be captured. The code is 
also able to perform calculations with large Lorentz fac- 
tors (r > 100), which allows us to probe astrophysically 
relevant regions of parameter space. 

2. METHOD 
We solve the relativistic hydrodynamics equations: 

^0=0 = (i) 

cy (p + 4P)uV + Pgr) = (2) 

where p is the proper mass density, P is the pressure, 
and u is the four velocity. We use units for which c = 1. 
We also evolve a passive scalar field, X, according to 

a ll (pXvf i ) = o. (3) 

This allows us to track the growth of the instability, de- 
marking which fluid elements arc in the ejecta and which 
are in the circumburst medium. We set X = for the 



ejecta behind the contact discontinuity, and X = 1 for 
the ambient gas in front. 

2.1. Numerics 

The numeric al method we employ is a v ariant of the 
TESS method (|Duffell fc ^ MacFadyen l20ll with specific 
application to problems involving rapid radial outflow. 

TESS is a finite- volume, moving-mesh hydrodynamics 
method which in its original form constructs its numer- 
ical mesh from a Voronoi tessellation of the computa- 
tional domain. The computational zones move with the 
local fluid velocity, and change their shape and size as 
they move and shear past one another. In this sense, the 
method is effectively Lagrangian. 

This idea can be improved if the moving mesh is not 
generated from a Voronoi tessellation, but instead is 
specifically tailored to a given flow. Most of the numeri- 
cal benefits from TESS do not come from the Voronoi tes- 
sellation, but from the simple fact that cells are allowed 
to move. This enables resolution of contact discontinu- 
ities to very high precision, and potentially allows for 
much longer timesteps. By moving the mesh, but choos- 
ing a particular shape for grid cells (other than Voronoi), 
the code can be adapted to the problem at hand. 

The TESS method defines its numerical mesh ab- 
stractly, so that all that is required to take into account 
a given mesh topology are the positions and volumes 
of computational zones, and the positions and areas of 
"faces" , which are simply defined to be any boundary be- 
tween two zones. This means more specific mesh topolo- 
gies can be generated without changing anything funda- 
mental about the method. 

We have previously made this adaptation for Keplerian 
flow in a gaseous disk (jDuffell fc MacFadven I [20121) . In 
that case, instead of Voronoi cells, the cells were chosen 
to be wcdgc-likc annular segments like those of a stan- 
dard polar grid, and they rotated with the local orbital 
velocity. We now specifically refer to the shearing-disk 
version as the DISCO code. 

In this work, we take a similar approach, but instead 
adapt the grid to a radial outflow, such as a jet. The 
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Fig. 2. — The numerical grid shears with the radial velocity of 
the gas. Each radial track moves independently, each behaving 
essentially as a ID lagrangian code. The upper panel shows the 
motion of grid cells, and the lower panel shows the shearing grid 
being employed to follow a jet-like explosion. 

wedge-like annular shape is still employed, but instead 
of rotating, the grid expands and shears radially (Figure 
[2]). The numerics are nearly identical to DISCO; we have 
merely changed the direction of grid motion. By moving 
the mesh radially instead of azimuthally, the code effec- 
tively behaves as a series of one-dimensional Lagrangian 
codes, coupled laterally by transverse fluxes. The mov- 
ing grid helps to capture the growth of the instability, 
and it also allows accurate calculations for large Lorentz 
factors T > 100. 

We integrate the equations in two dimensions assuming 
axisymmctry. The computational domain consists of an 
angular wedge < 9 < w/16. We use an angular resolu- 
tion Ng — 1600 and a resolution in the radial dimension 
of Nr ~ 12800, which were chosen for a resolution of 
A9 ~ Ar/r ~ 10~ 4 . The number of radial zones varies 
throughout the course of the calculation, because zones 
are added and removed from the computational domain 
as necessary during the course of the calculation. The 
aspect ratio of a typical computational zone is of order 
unity, but if a cell becomes too long or short, the grid is 
dynamically refined or de-refined by splitting or merging 
zones. Aspect ratios in general range between 5 and 1/5. 
The inner and outer radial boundaries of the domain can 
also move during the time integration, so that only frac- 
tion of the dynamical range must be covered at any one 
time. 
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Fig. 3. — The power spectrum of the Raylcigh- Taylor for the 
three examples shown in Figure [T] Curves are calculated by mea- 
suring the Fourier transform of 8r/r in front of the Raylcigh- Taylor 
instability. There is a break in the power spectrum at the causality 
scale A0 ~ 1/r. 

2.2. Initial Conditions 

The calculation begins with a hot explosion with a 
given energy E and mass M with characteristic Lorentz 
factor r = E/M, within a negligibly small radius, in a 
cold uniform-density medium with density pq. This spec- 
ifies the problem completely, provided we have a defini- 
tion for "negligibly small radius" . The initial radius Rq 
of the explosion must be small enough that the "coast- 
ing phase" at t ~ TR and the "spreading phase" at 
t ~ T 2 i?o occur long before the shell begins to deceler- 
ate, at the deceleration time 



In this limit (T 2 Rq << i 7 ), the initial structure of the 
fireball has no effect on what happens after t > £ 7 (we 
have checked this empirically). This is the so-called 
"thin-shell" limit. At t = t 7 , the shell decelerates, which 
causes the formation of the forward and reverse shocks. 
This generates a self-similar solution which persists until 
t = tM, where 

t M = (M/po) 1/3 - (5) 

During this phase of expansion, between i 7 and tM, the 
Ray leigh- Taylor instability is most significant. 

To capture the first phases of evolution, we evolve the 
explosion in one dimension from the initial fireball until 
a time t < t 7 , when the mass in the forward shock is 
still a small fraction of the mass of the ejecta. Then the 
code is switched to evolve the system in 2D, to capture 
the growth of the instability. At this point, we also set 
the value of X on either side of the blast wave (X = 
behind and X = 1 in front). The initial parameters are 
E, M, and po, but due to scale invariance, in this thin- 
shell limit the only parameter which must be varied is 
the energy per unit mass in the fireball, E/M = F. Thus 
all the parameter space of this problem is reduced to this 
single variable. In this study, we use values of T = 10, 
30, and 100. We have also performed calculations with 
F = 300 but resolution was insufficient to correctly evolve 
the instability. The T > 100 case will be presented in a 
future work. 

In principle, the instability grows from some seed per- 
turbations which break the spherical symmetry of the 
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Fig. 4. — The same calculation has been performed in 2D assum- 
ing axisymmetry, and in ID assuming spherical symmetry. The 
instability is present in two dimensions, and has a clear impact on 
the spherically averaged fields. This is a snapshot of proper density 
at t = tM for r = 30, comparing the ID and 2D solutions. The 
top panel shows the scatter of density present in two dimensions, 
compared to the two-shock solution with sharp contact discontinu- 
ity found in ID. The second panel shows an angular average of the 
two-dimensional field, which demonstrates that the contact discon- 
tinuity is completely disrupted. Panels 3-5 show a logscale profile 
in density, pressure and Lorentz factor, showing clearly the posi- 
tion of the reverse shock, which is pushed back by the insta bility . 
The bottom panel shows the turbulent energy fraction, ert ( 33.20 . 



initial conditions. However, we found that for weak 
enough perturbations the saturated instability does not 
depend (in a statistical sense) on how the perturbations 
are seeded. This includes the limit where there are no 
explicit perturbations, but the instability is seeded with 
numerical noise, which is extremely small but still large 
enough to break the symmetry. All calculations per- 
formed here did not use explicit perturbations to seed 
the instability. 
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Fig. 5. — Another comparison between the one-dimensional and 
two-dimensional versions of this experiment, for Y = 30. In this fig- 
ure we plot the positions of the characteristic waves present within 
the shell. The forward shock is unaffected by the instability, the 
contact discontinuity is sheared out into the large shaded region, 
and the reverse shock is pushed away from the forward shock. 

3. RESULTS 

In Figure [T] we demonstrate how relativistic effects 
can affect the character of Rayleigh- Taylor turbulence. 
We show a snapshot at t = t^j for various values of 
r = E/M, which is roughly a proxy for the initial Lorentz 
factor. For larger Lorentz factors, the fastest growing 
modes can be found on smaller angular scales. This 
has to do with causality; the saturated instability favors 
large-scale modes, but there is a large-scale cutoff given 
by the horizon scale, A9 ~ 1/T. Modes larger than this 
scale do not grow. This is most explicitly shown in the 
power spectrum of the contact discontinuity, which we 
have plotted in Figure [3j The perturbation also propa- 
gates further both upstream and downstream for larger 
values of T. For the turbulence moving back into the 
ejecta this is due to increased pressure and density gradi- 
ents in the unstable region. For the forward-propagating 
turbulence, it is due to the fact that the gas is compressed 
into a thin shell ~ r/T 2 behind the forward shock, and 
the width of this shell decreases with T. 

3.1. ID vs 2D 

The primary effects of the Rayleigh- Taylor instability 
on the shock dynamics are displayed in Figure [4] This is 
a snapshot of the Shockwave at t = tM for T = 30. We 
show the radial profile and compare it with the same cal- 
culation performed in ID assuming spherical symmetry. 
In the ID case, we see the standard two-shock solution 
which has been extensively used to describe GRB phe- 
nomenology 

From this figure, the most obvious distinction between 
one and two dimensions is the complete disruption of the 
contact discontinuity in 2D, which of course is expected 
since it is the contact discontinuity which is unstable. 
The second plainly visible distinction here is in the re- 
verse shock. Clearly in the 2D case, the instability has 
managed to perturb the reverse shock, and this has al- 
tered the shock's characteristic speed. In the 2D case, the 
reverse shock propagates much further back than in ID. 
This apparently does not happen to the forward shock, as 
it remains undisturbed. It is possible that this conclusion 
could be challenged by a three dimensional calculation, 
but we do not expect this to be the case, as the instability 
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Fig. 6. — The fireball converts itself from a thin shell into the two-shock structure, and as t — > tM the characteristic waves obey 
self-similarity. This remains true for the unstable 2D version of the problem. Here we have plotted the position of the characteristic waves 
in terms of the self-similarity variable, \ (Eq. [fij. Shortly after t = tp,j, the characteristic waves fall behind the shock. In the upper right 
panel we show the fraction of energy which is perturbed by the instability, roughly 10% at t = t]\f. In the lower right panel we show how 
the position of characteristic waves scales with Lorentz factor. 



would dissipate more efficiently in 3D. 

Characteristic waves are shown in Figure [5l where we 
compare ID and 2D shock dynamics. Here again, we 
can see important effects of the instability: What was 
a contact discontinuity becomes smeared over a larger 
region in 2D, overtaking the reverse shock but leaving 
the forward shock alone. 

We can compare different values of T on a similar foot- 
ing by expressing the positions of characteristic waves in 
terms of the self-similarity parameter, 



X 



1 - r /t 

1 f shock I ^ 



(6) 



We plot the value of x f° r a U characteristic waves in 
Figure E] for T = 10,30, and 100. As t -> t M , the char- 
acteristic waves approach constant values of \- From 
these three examples, it appears that the forward posi- 
tion of the Rayleigh- Taylor fingers approaches a value of 
X ~ 2.5 which is independent of T. We might expect 
from this that turbulence should be present in about one 
third of the energy behind the shock, though in prac- 
tice we find this fraction is closer to 10%. This is the 
percentage which is turbulent and therefore magnetized 
by Rayleigh- Taylor (we calculate the magnetization in 
^3.2|) . After t = tu, the shock moves away from the fin- 
gers quickly and have no further impact on the forward 
shock dynamics. 

3.2. Turbulence 

It will be useful to measure the kinetic energy of large 
eddies in the turbulent region of the solution, since that 
will provide insight into the turbulent cascade which 
would take place in a highly-resolved 3D calculation. It 
turns out that the relevant quantity, the turbulent energy 
fraction, 



i-RT 



Turbulent Kinetic Energy 
Thermal Energy 



(7) 



can be measured in a relatively simple manner. In our 
analysis, we have converted 2D data into ID data by 
averaging in a conservative way. That is, for a given ra- 
dial annulus, the total mass, energy, and momentum of 
that annulus are calculated, and from these quantities 



the conservatively averaged density, pressure, and four 
velocity are determined (this is how the middle panels 
of Figure H] were generated). We note that this averag- 
ing process typically results in a slightly larger pressure 
than would result from arithmetically averaging the pres- 
sure by volume on an annulus. This is because when the 
momentum of an annulus is calculated, there is cancella- 
tion due to the variations in velocity with angle, and this 
"lost" kinetic energy is interpreted as thermal energy in 
the averaged quantities. 

In other words, the excess of thermal energy from av- 
eraging conservatively over angle is exactly the turbulent 
kinetic energy contained in that annulus. This provides 
a very simple means of calculating the fraction of energy 
in large eddies: 

tRT = (8) 

where AP is the difference between "conservatively av- 
eraged" pressure and volume-averaged pressure. 

We have plotted a snapshot of e^T in the lower panel 
of Figure |U There is clearly variation with radius, but in 
this snapshot it roughly peaks around e^T ~ 0.02. For 
the three cases studied, this turbulent fraction peaks at 
time t ~ tM with this same peak value, roughly indepen- 
dent of r. 

This turbulent energy is free to cascade to smaller 
scales (as it would in a 3D calculation) and to amplify 
magnetic fields (as it would if magnetic fields were in- 
cluded in our calculation). In fact, it is straightfor- 
ward to estimate the saturated ma gnetic energy frac 
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tion, because it has been shown ([Zhang et al 
iZrake fc MacFadven |[20T3) that it is roughly in equipar- 
tition with the turbulent kinetic energy: 

e B ~ CRT ~ 10- 2 . (9) 

This is the magnetization induced by Raylcigh- Taylor 
alone, though magnetic fields could be amplified by other 
mechanisms. 

3.3. Afterglow Emission 

We have calculated light curves in one and two di- 
mensions, in order to explore the observational conse- 
quences of the instability. Relativistic Doppler effects 
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Fig. 7. — The reverse shock is directly observed in the early 
afterglow emission. We plot microwave light curves for T = 30 for 
the ID (stable) and 2D (unstable) versions of the problem. The 
second peak is due to the forward shock, which is unaffected by 
the instability. The first peak is emission from the reverse shock. 
This flash occurs at a later time and the peak is narrower. 

are responsible for the hard gamma-ray and X-ray af- 
terglow emission, but the reverse shock generally has a 
much lower Lorentz factor than the forward shock, so it 
cannot compete with the forward shock at high frequen- 
cies. At lower frequencies (optical, microwave, radio) the 
flux from the forward and reverse shocks can be compa- 
rable, so in these wavelengths we search for observational 
signatures of the Rayleigh- Taylor instability. 

By implementing a simple model for synchrotron emis- 
sion directly into the hydrodynamical evolution, we cal- 
culate a microwave light curve for our T = 30 model. We 
calculate the flux from a given computational zone at a 
given observer time from the formula: 



F oc 



P B 



7 2 (1 



;Q0>) 



(10) 



where the magnetic energy is assumed to be proportional 
to the thermal energy, 

BocVP (11) 
and frequency dependence assumes "slow cooling" 



l(l"P)/2 



V < V r , 

V > V„ 



(12) 



where we choose p = 2.5, and where v m is the character- 
istic synchrotron frequency, 



B{P/pf 



(13) 



This is a simplif i ed ve rsion of the emission model in 
iVan Eerten etaTl (|20Tot) . 

We show the light curve in Figure [7] for T = 30. To 
single out signatures of Rayleigh- Taylor, we compare ID 
and 2D results. The forward shock emission is unaf- 
fected, as the forward shock is not perturbed by the in- 
stability. The reverse shock emission is affected by the 
instability; the peak magnitude is not diminished signif- 
icantly, but the peak occurs at a later time in 2D than 
in ID. 

4. SUMMARY 

We have computed the development and saturation of 
Rayleigh- Taylor instability in relativistic fireballs using 
a moving computational mesh. The instability does not 
perturb the forward shock, but it completely disrupts 



the contact discontinuity between the ejecta and the sur- 
rounding medium, and it significantly modifies the dy- 
namics of the reverse shock, which may cause reverse 
shock emission to peak at a later time. The relativistic 
shear flow generated by Rayleigh- Taylor fingers may also 
provide a site for cosmic ray acceleration. 

The Rayleigh- Taylor instability is effective in stirring 
up turbulence; we find that turbulent kinetic energy 
reaches about 1% of thermal energy wherever the tur- 
bulence is present, which should rapidly amplify mag- 
netic fields up to 1 % of thermal by the oper ation of 
small-scale dynamo (jZrake fc MacF adven II2013T) during 
the deceleration phase (between i 7 and 4m)- The turbu- 
lence is present in a region behind the forward shock, but 
Rayleigh- Taylor fingers attempt to reach out towards this 
shock, penetrating about 10% of its energy. The reverse 
shock is impacted immediately, and it is pushed away 
from the forward shock, strengthening it and causing its 
emission to be observed at a later time. This motivates 
3D calculations, as these numbers might change when 3D 
turbulence effects are taken into account. 

This emphasizes the importance of accounting for fluid 
instabilities when considering the dynamics and obser- 
vational signatures of gamma ray bursts. Figure |3] gives 
some hope that an effective ID description of the dynam- 
ics may be possible; adding some kind of diffusive term 
in ID to shocked fluid elements in the ejecta, for exam- 
ple, might reproduce physical profiles close to those in 
Figure |3J 

This motivates a more comprehensive study which in- 
cludes more general circumburst profiles, larger Lorentz 
factors (converged results with T ~ 300 or higher are 
much more computationally intensive but certainly pos- 
sible), and potentially 3D turbulence effects. Because 
of scale invariancc, our only free parameter is T, the 
characteristic Lorentz factor. The dynamics and syn- 
chrotron emission both can be re-scaled to arbitrary ex- 
plosion energy and circumburst density, which means 
that a suite of calculations for different values of V will 
compl etely cover the parameter space o f all such explo - 
sions ([Van Eerten fc MacFadven I I2012t IGranot 1 120121) . 
This motivates us to implement a detail ed synchrotron 
emission model (|Van Eerten et ai"ll2010[) . to perform a 
comprehensive study of emission from GRB outflows. All 
of this is planned as future work. 
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